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Abstract:  In  this  study  we  describe  a  numerical  algorithm  which  deduces 
characteristics  of  the  bottom  bathymetry  given  free  surface  elevation  records  dense  in 
space  but  sparse  in  time.  The  method  makes  use  of  the  Levenberg-Marquardt  numerical 
optimization  scheme  in  conjunction  with  a  time-domain  nonlinear  model.  Iteration 
occurs  until  the  mismatch  between  the  free  surfaces  of  the  data  and  model  are 
minimized;  the  bathymetry  is  adjusted  in  order  to  achieve  this  minimum.  The 
sensitivity  measure  is  a  by-product  of  the  calculation,  and  determines  the  invertibility 
of  the  system.  Due  to  robustness  concerns,  we  limit  ourselves  to  deduction  of 
bathymetric  profile  parameters.  Tests  of  the  system  using  monochromatic,  irregular  and 
groupy  waves  show  favorable  results;  the  latter  is  particularly  notable  give  the 
difficulty  standard  inversion  methods  have  had  with  groupy  waves.  A  two-stage  system 
is  also  outlined,  in  which  a  simple  parameterization  for  a  nearshore  bar  is  developed 
and  utilized.  The  first  stage  determines  the  mean  profile,  while  the  second  stage 
determines  the  bar  characteristics  using  the  first  stage  results  as  the  initial  iterate.  To 
extend  the  method’s  capabilities  further,  the  use  of  phase  speed  records  is  discussed. 


Introduction 

Nearshore  wave  and  circulation  models  have  become  much  more  sophisticated 
in  both  structure  and  physics.  However,  the  accuracy  of  these  models  is  highly 
dependent  on  the  quality  and  resolution  of  the  nearshore  bathymetry.  Direct 
measurement  of  the  bathymetry  is  difficult  under  the  most  benign  conditions. 
Additionally,  the  nearshore  waves  and  currents  obviously  react  to  the  changes  in  the 
bathymetry  in  deterministic  ways.  The  physical  relationships  between  free  surface 
evolution  and  bathymetric  variations  (represented  in  numerical  models  of  wave 
evolution)  can  be  exploited  to  deduce  the  bathymetry  given  realizations  of  the  free 
surface  or  other  properties  (e.g.,  phase  speeds,  etc.). 

Concurrently,  remote  sensing  capabilities  have  also  undergone  significant 
development.  Synthetic  aperture  radar  (SAR),  video,  electro-optical  sensors,  and  X- 
band  (marine)  radar  are  capable  of  high  resolution  measurements  of  light  intensity  from 
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the  free  surface.  Wave  kinematics  can  thus  be  calculated  by  linking  variations  in 
intensity  to  phase  lines  of  the  waves.  Direct  sensing  of  the  free  surface  or  surface 
currents  can  be  accomplished  using  airborne  lidar  (for  a  single  snapshot  of  the  free 
surface)  or  interferometric  synthetic  aperture  radar  (INSAR).  These  latter  systems  will 
give  either  a  direct  free  surface  measurement  (lidar)  or  an  inferred  measurement  from 
the  currents  (INSAR).  The  challenge,  therefore,  lies  in  incorporating  these  data  into 
bathymetry  deduction  algorithms. 

The  problem  of  bathymetry  deduction  from  remotely-sensed  wave  characteristics  is  one 
that  has  been  a  military  focus  since  World  War  II,  but  has  received  much  attention 
recently.  These  recent  developments  can  be  roughly  classified  into  two  areas.  The  first 
approach  focuses  on  gleaning  wave  kinematics  and  dynamics  from  the  free  surface 
imagery,  then  using  the  linear  dispersion  relation: 

co  =  gk  tanh  kh  (1) 

in  which  co  is  the  radian  wave  frequency,  g  is  gravitational  acceleration,  and  k  is  the 
wavenumber  (=2jdL,  where  L  is  the  wavelength).  Bell  (1999)  used  phase  speeds 
calculated  from  X-band  radar  imagery  and  Equation  (1)  to  infer  the  bathymetry,  with 
favorable  results.  Holland  (2001)  performed  a  similar  analysis  using  phase  speeds  from 
video  data,  while  Dugan  et  al  (2001)  utilized  kinematics  from  airborne  electro-optical 
systems. 

The  second  approach  combines  data  and  numerical  models  into  a  numerical 
optimization  problem,  in  which  the  total  mismatch  between  the  data  and  the  model 
results  is  iteratively  minimized,  with  the  bathymetry  correspondingly  adjusted.  Kennedy 
et  al.  (2000)  and  Narayanan  and  Kaihatu  (2000)  are  examples  of  this  approach.  We  also 
note  that  Grilli  (1997)  created  nomographs  for  direct  calculation  of  the  bathymetry  from 
phase  speeds  based  on  results  from  his  fully  nonlinear  boundary  element  model. 

Our  development  uses  the  second  approach,  and  is  a  refinement  of  that  detailed 
in  Narayanan  and  Kaihatu  (2000)  in  that  we  are  now  able  to  “zoom  in”  on  nearshore 
features  such  as  sandbars.  Additionally,  we  explicitly  investigate  the  improvement  of 
this  model  over  direct  inversion  of  Equation  (1)  with  regards  to  wave  groups.  We 
conclude  with  considerations  of  future  directions  for  this  work. 

Numerical  Models  and  Data 

In  the  interest  of  computational  efficiency  during  the  initial  development  of  the 
bathymetric  deduction  algorithm,  we  used  the  Korteweg-deVries  (KdV)  equation 
(Korteweg  and  deVries  1895)  as  the  wave  model.  Throughout  this  study,  we  will  be 
using  output  from  the  KdV  model  as  synthetic  data  for  system  testing;  any  physical 
shortcomings  contained  in  the  wave  model  is  therefore  reflected  in  the  data  as  well,  and 
so  thus  is  consistent.  The  modular  nature  of  the  completed  algorithm  is  such  that  any 
wave  model  can  be  used  instead  of  the  KdV  model. 

We  assume  that  the  free  surface  imagery  is  in  the  fonn  of  “snapshots”  of  the  free 
surface,  dense  in  space  but  well-separated  in  time.  This  sort  of  free  surface  data  stream 
would  be  akin  to  having  a  lidar-equipped  aircraft  make  several  sweeps  over  the  same 
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area.  This  is  in  fact  a  severe  test  of  the  algorithm,  since  the  time  span  between  snapshots 
is  greater  than  a  dominant  wave  period.  The  likelihood  is  great  that  a  more  generous 
data  set  will  yield  better  results  for  the  bathymetry. 

Bathymetric  Parameterizations  and  Optimization  Techniques 

The  numerical  optimization  technique  used  herein  is  the  Levenberg-Marquardt 
method,  which  is  an  overrelaxed  nonlinear  least  squares  technique  (see  Press  et  al.  1986 
for  sample  algorithms).  This  method  has  been  used  in  many  applications,  and  excels  in 
the  case  of  highly  nonlinear  optimization  with  multiple  local  minima  and  elongated, 
shallow  error  surfaces. 

At  the  time  of  writing,  the  authors  are  preparing  a  journal  manuscript  on  the 
topic  at  hand;  the  full  technique  is  explained  in  an  appendix  of  the  manuscript.  In  the 
interest  of  brevity,  we  simply  write  the  matrix  equation  to  be  solved: 

Ah  =  (ATA  +  ^2I)ATAr/  (2) 


where  Ah  is  a  vector  containing  the  difference  between  the  prior  and  new  iterates  of 
water  depth;  I  is  an  identity  matrix,  A rj  is  a  matrix  of  the  data -model  mismatches,  //  is 
an  over-relaxation  parameter,  and  the  matrix  A  is  known  as  the  “sensitivity”  matrix: 
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given  m  total  parameters  (water  depths  to  be  determined,  in  this  case)  and  n  total  data 
points.  These  data  could  be  distributed  in  space  (as  is  the  case  for  wave  staffs  measuring 
time  series  at  discrete  locations)  or  time  (remotely  sensed  snapshots  of  the  free  surface 
taken  at  discrete  times).  Multiplying  by  the  transpose  creates  an  mxm  square  matrix. 
The  diagonal  elements  of  this  square  matrix  comprise  an  mx  1  vector: 
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Each  element  of  this  vector  is  known  as  a  sensitivity  measure,  and  is  a  key  to  the 
solvability  of  the  system.  In  general,  the  relative  amplitudes  of  all  the  m  sensitivity 
measures  must  be  within  three  orders  of  magnitude  in  order  for  the  system  to  be 
invertible  and  the  parameters  identifiable.  This  insures  that  the  search  for  any  one 
parameter  does  not  dominate  the  iteration  process.  The  sensitivity  measures  are  not  the 
only  considerations  present  in  the  inversion  process.  However,  the  sensitivity  measures 
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are  also  useful  in  other  aspects  of  the  inversion  problem,  such  as  array  design  and  the 
formulation  of  data  sampling  strategies.  We  will  discuss  this  aspect  further  in  a  later 
section. 

During  the  initial  development  of  this  system,  we  endeavored  to  be  able  to  find  an 
arbitrary  number  of  water  depths  in  our  domain,  at  least  as  many  as  we  had  grid  points 
in  the  wave  model.  However,  difficulties  with  nonconvergence  caused  us  to  consider 
using  various  parameterized  forms  of  the  bathymetric  cross-shore  profile,  with  the 
various  free  parameters  within  the  assumed  forms  (rather  than  the  depths  themselves) 
becoming  the  targets  of  the  search.  We  considered  several  parameterized  forms, 
including  Dean’s  equilibrium  beach  profile  (Dean  1977),  exponential  forms  (Bodge 
1992),  and  polynomial  forms.  We  detennined  that  both  Dean’s  equilibrium  beach 
profile  and  exponential  forms  worked  best  with  the  inversion  scheme  detailed  here.  In  a 
later  section  we  describe  how  the  ostensibly-limiting  form  for  bathymetry  can  be  used 
as  an  initial  guess  for  more  detailed  inversions  in  the  nearshore. 

Monochromatic  and  Random  Waves 

Our  tests  of  the  system  with  both  monochromatic  and  irregular  waves  have  been 
detailed  in  Narayanan  and  Kaihatu  (2000).  We  used  Dean’s  equilibrium  beach  profile 
form  h(x)=Axm  ,  and  use  the  method  described  in  the  previous  section  to  determine  the 
free  parameters  A  and  m .  A  set  of  free  surface  snapshots  would  be  generated  by 
propagating  the  waves  over  bathymetry  for  which  the  parameters  are  known,  then 
sampling  the  modeled  free  surface  elevations  over  the  entire  domain  at  three  discrete 
points  in  time  (t=  16,  20  and  24  seconds).  The  domain  was  35  meters  long,  and  the  wave 
had  a  peak  period  of  5  seconds.  The  bathymetric  parameters  were  then  varied  and  the 
iteration  procedure  begun;  each  step  of  the  iteration  yielded  trial  values  for  the 
bathymetric  parameters  A  and  in .  Convergence  of  the  system  to  the  actual  parameter 
values  was  highly  dependent  on  initial  guess.  Plots  of  the  results  for  the  random  wave 
case  can  be  seen  in  Narayanan  and  Kaihatu  (2000).  An  empirically-derived  map  of 
convergence  as  a  function  of  initial  position  in  parameter  space  is  shown  in  Figure  1. 
The  plot  reveals  that  initial  guesses  which  tend  to  overestimate  the  offshore  depth  (high 
values  of  A)  tend  to  lead  to  convergence  more  than  those  parameter  values  which 
underestimate  the  offshore  depth. 

Wave  Groups 

In  our  investigations  into  the  applicability  of  direct  inversion  methods  utilizing 
the  linear  dispersion  relation  (Equation  1),  it  was  found  that  the  case  of  wave  groups 
required  special  treatment.  Both  monochromatic  and  random  wave  trains  have  free 
surface  undulations  which  are  more  or  less  strongly  coupled  to  variations  in  the 
bathymetry.  Strong  wave  groups  (those  comprised  of  a  small  number  of  frequency 
components)  have  free  surface  characteristics  which  may  be  more  due  to  the 
interference  patterns  inherent  in  the  group  rather  than  bathymetric  variations.  This 
feature  tends  to  cause  underestimation  of  wavenumbers  and  results  in  overprediction  of 
the  water  depth. 
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Figure  1.  Empirical  convergence  map,  Dean’s  equilibrium  beach  profile  under 
random  waves.  Initial  guesses  of  A  and  m  corresponding  to  open  circles  will 
diverge,  while  those  corresponding  to  asterisks  will  converge. 

Consideration  of  additional  data  in  time  could  be  expected  to  improve  results. 
This  is  true  of  random  waves,  since  the  peak  spectral  values  gleaned  from  the  data  and 
used  in  the  linear  dispersion  relation  for  inversion  would  be  more  statistically 
representative  of  the  actual  wave  field.  However,  this  does  not  apply  to  strong  wave 
groups.  Figure  2  shows  the  results  of  a  numerical  experiment  detailed  in  Kaihatu  et  al. 
(1999),  in  which  a  wave  model  was  used  to  propagate  random  and  groupy  waves  over  a 
bathymetric  cross-shore  transect  measured  at  the  US  Army  Field  Research  Facility  at 
Duck,  NC.  Over  2000  temporal  realizations  of  the  free  surface  (amounting  to  around 
200  seconds  of  time)  were  generated.  A  wavelet  transform  was  used  to  glean 
wavenumbers  for  each  point  along  the  cross-shore  transect  for  each  temporal 
realization,  then  averaged  across  all  realizations.  This  gave  a  statistically  significant 
peak  wavenumber  for  each  point  along  the  cross-shore  transect.  This,  along  with  the 
peak  frequency  from  the  incident  spectrum,  was  used  with  the  linear  dispersion  relation 
to  solve  for  the  water  depth.  Figure  2a  shows  the  case  for  random  waves;  despite  the 
crudity  of  the  inversion  method,  the  depths  match  quite  well.  Figure  2b  shows  the  result 
of  the  same  analysis  using  wave  groups.  Despite  the  significant  temporal  averaging,  the 
strong  group  signature  is  still  apparent  in  the  result. 

The  inversion  procedure  described  earlier  was  then  applied  to  the  case  of  wave 
groups.  Results  are  shown  in  Figure  3.  It  is  apparent  that  the  procedure  works  quite 
well.  Figure  3a  shows  the  initial  and  true  profiles  under  wave  groups,  while  Figure  3b 
shows  the  empirical  convergence  map.  There  appears  to  be  a  cluster  of  convergence 
near  the  actual  value  of  A  and  in  within  which  initial  guesses  for  A  and  m  will  lead  to 
convergence.  The  success  of  the  procedure  for  wave  groups  indicates  that  the  benefit  of 
this  method  is  the  lack  of  wave  kinematics  or  other  assumptive  properties  that  require 
extraction  from  imagery. 
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Figure  2.  Results  of  direct  depth  inversion  using  linear  dispersion  relations  and 
wave  kinematics  gleaned  from  numerical  results  by  wavelet  transforms.  Solid  line: 
estimated  bathymetry  from  analysis.  Dashed  line:  actual  bathymetry  from  Duck, 
NC.  Left  (a):  Random  waves.  Right  (b):  Wave  groups. 
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Figure  3.  Results  of  bathymetric  inversion  using  wave  groups.  Left  (a):  Initial  and 
final  bathymetric  profiles.  Right  (b):  Empirical  convergence  map 

Two-Stage  Inversion 

We  further  refined  the  bathymetric  inversion  scheme  by  adding  a  second  stage 
to  the  iteration  routine.  This  second  stage  uses  the  estimation  of  the  mean  profile  as  the 
first  guess  for  a  second  round  of  iterations. 

As  with  the  single  stage  inversions  described  earlier,  we  use  a  parameterized 
bathymetry  for  the  second  stage  (as  with  the  previously-described  inversion,  we 
encountered  convergence  problems  when  we  attempted  to  use  every  water  depth  in  our 
numerical  grid  as  a  parameter  for  which  to  search).  Our  parameterization  of  this  second 
stage  bathymetry  includes  a  simplistic  formulation  to  account  for  the  bar: 
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h2  (x)  =  a  sin  J3(x  -  x0  ) 


(5) 


where  the  subscript  “2”  denotes  the  results  of  the  second  stage  searchm  a  and  ft  are 
shape  factors  for  the  bar,  and  xo  is  the  location  of  the  offshore  edge  of  the  bar.  In  theory 
one  could  search  for  all  five  parameters  in  one  stage;  however,  a  sensitivity  analysis  had 
determined  that  the  search  for  the  bar  parameters  must  be  perfonned  over  a  much 
smaller  section  of  the  domain  than  that  of  the  mean  profile  parameter  search.  Thus,  the 
search  takes  place  in  two  sets  of  iterations.  The  first  set  of  iterations  yields  the  mean 
profile  parameters  A  and  m,  and  the  second  search  seeks  the  other  three  parameters. 

We  tested  the  two-stage  inversion  system  with  bathymetry  at  Duck,  NC.  This 
particular  set  of  bathymetry  was  collected  during  the  October  1990  DELILAH 
experiment,  during  which  a  bar  was  almost  always  present.  A  TMA-type  wave 
spectrum  was  run  over  this  bathymetry  and  the  free  surface  information  along  a  cross¬ 
shore  transect  of  the  domain  stored  at  t=  16,  20  and  24  seconds.  This  data  then  seeded 
the  inversion.  Ligure  4  shows  the  result  of  the  two-stage  inversion.  The  first  stage 
parameters  do  well  in  replicating  the  mean  profile  at  Duck  for  much  of  the  domain,  and 
even  capture  the  mean  of  the  bar  fairly  well.  The  second  stage  results  (shown  in  the 
extreme  nearshore  area  of  the  domain)  captures  some  of  the  reverse  curvature  of  the  bar. 
As  mentioned  before,  the  sensitivity  analysis  for  this  scenario  indicated  that  the  offshore 
extent  of  the  second  stage  iterations  needed  to  be  very  close  to  shore  if  the  inversion 
was  to  be  sufficiently  robust.  This  makes  intuitive  sense,  since  one  would  not  expect 
information  well  offshore  of  the  bar  to  be  essential  to  the  detennination  of  bar 
parameters;  however,  it  is  reassuring  that  a  quantitative  sensitivity  analysis  confirms  a 
qualitative  observation. 


Figure  4.  Results  of  two-stage  inversion  process:  comparison  to  Duck  bathymetry. 
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Phase  speed 

Though  the  overall  approach  described  herein  is  powerful  and  has  wide  potential 
applications,  there  are  several  obvious  disadvantages  with  our  specific  application.  The 
reliance  on  bathymetric  parameterization  is  a  severe  limitation,  and  would  hamper  an 
extension  to  longshore  variations  of  bathymetry.  Additionally,  the  use  of  temporally 
sparse,  spatially  dense  free  surface  elevations  as  data  limits  potential  applicability  of  the 
system  to  a  small  range  of  remote  sensing  instrumentation  (interefereometric  SAR, 
lidar).  Adaptation  of  the  inversion  algorithm  to  phase  speed  measurements  would 
greatly  expand  the  possible  remote  sensing  platforms  which  can  be  used  in  tandem,  as 
mentioned  earlier. 

It  can  be  somewhat  cumbersome  to  calculate  phase  speeds  from  a  free  surface 
time  domain  wave  model,  though  correlation  algorithms  such  as  those  used  by  Bell 
(1999)  and  Holland  et  al.  (2001)  can  be  used.  However,  it  is  more  expedient  to  use  a 
frequency-domain  model  for  this  application,  since  it  is  simpler  to  calculate  phase 
speed,  even  with  nonlinear  effects. 

To  investigate  the  efficacy  of  calculating  phase  speed  from  a  frequency  domain 
model,  we  use  the  model  of  Kaihatu  and  Kirby  (1995),  with  extensions  by  Kaihatu 
(2001),  to  evolve  a  wave  spectrum  over  a  planar  bathymetry.  This  model  is  developed 
in  terms  of  the  complex  amplitude  of  the  free  surface,  with  resonant  triad  interactions 
representing  the  shallow  water  wave  nonlinearity.  The  particular  advantages  of  this 
model  are  addressed  in  Kaihatu  and  Kirby  (1995).  Substituting  a  real-valued  amplitude 
and  phase  in  place  of  the  complex  amplitude,  then  retaining  the  imaginary  part  of  the 
resulting  equation,  yields  a  direct  expression  for  the  phase  speed.  We  initialize  the 
model  with  a  TMA-type  spectrum  (peak  period  of  X  seconds  and  a  gamma  of  Y),  and 
run  the  spectrum  over  a  sloping  bathymetry  with  an  offshore  depth  of  8  meters  and  a 
bottom  slope  of  0.0005.  Figure  5  shows  the  nonlinear  phase  speed  calculation  from  the 
model  (using  the  peak  frequency),  with  a  corresponding  calculation  for  linear  phase 
speed.  The  nonlinear  phase  speed  is  greater  than  that  from  linear  theory  for  the  offshore 
part  of  the  domain.  As  breaking  begins  to  dominate,  the  waveheight  at  the  peak 
frequency  becomes  smaller  and  nonlinear  effects  are  reduced.  If  video  data  is  utilized, 
the  temporal  statistics  (due  to  the  long  measurement  time)  can  approach  those  of  the 
model,  thus  allowing  a  one-to-one  correspondence  (in  theory)  between  average 
characteristics  from  the  model  and  those  from  the  data.  Thus  phase  speed  estimates 
would  be  representative  of  the  same  amount  of  time,  lending  confidence  to  our  use  of 
phase  speed  for  depth  inversion. 

Summary 

We  have  outlined  a  method  for  determining  water  depth  characteristics  from 
remotely  sensed  data.  In  contrast  to  earlier  data-intensive  methods,  we  use  numerical 
optimization  and  a  numerical  wave  model  to  obtain  water  depth  estimates.  The 
numerical  technique  used  (Levenberg-Marquardt  optimization)  is  well-suited  to  the 
potential  nonuniqueness  of  the  inversion  process. 

Difficulties  with  convergence  compelled  us  to  utilize  bottom  profile 
parameterization  as  a  means  for  obtaining  bathymetric  characteristics.  We  obtained 
excellent  results  using  free  surface  elevations  gleaned  from  the  wave  model  as  “data,” 
and  using  the  numerical  algorithms  to  iterate  on  the  bathymetric  parameters.  This  was 
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true  for  monochromatic,  random,  and  groupy  waves.  The  latter  case  is  notable  in  that 
straightforward  depth  inversion  using  the  linear  dispersion  relation  (Equation  1)  proved 
to  be  problematic  since  the  free  surface  characteristics  are  more  a  function  of  the 
interference  patterns  of  the  waves  in  the  group  than  the  bathymetry  variations. 


Figure  5.  Phase  speed  estimates  from  the  nonlinear  frequency  domain  model  of 
Kaihatu  and  Kirby  (1995).  Phase  speed  calculated  from  the  peak  frequency  of  a 
wave  spectrum  shoaled  over  a  sloping  bottom.  Solid  line:  linear  phase  speed. 
Dashed  line:  nonlinear  phase  speed. 

We  then  also  investigated  the  parameterization  of  a  nearshore  bar  in  order  to 
extend  the  original  algorithm  by  implementing  a  two-stage  process.  The  mean  profile 
parameters  would  be  obtained  from  the  first  stage,  and  then  used  as  an  initial  iterate  for 
the  second  stage,  which  would  yield  the  bar  parameters.  We  tested  this  two-stage 
system  using  a  bathymetric  transect  from  Duck,  NC.  Both  the  mean  profile  and  the 
nearshore  bar  parameters  matched  the  actual  profile  fairly  closely. 

In  the  algorithm  testing  we  assumed  that  our  data  consisted  of  spatially  dense 
snapshots  of  the  free  surface  elevations  at  discrete  points  in  time.  This  focus  limits  the 
applicability  of  the  system  to  a  small  set  of  remote  sensing  instrumentation.  In  order  to 
overcome  this  limitation,  we  intend  to  focus  on  the  use  of  phase  speeds  from  both  the 
model  and  data  to  perform  the  inversion.  Phase  speed  estimates  are  obtainable  from 
many  more  remote  sensing  platforms  than  free  surface  elevations,  and  also  exhibit  less 
variability  than  free  surface  elevations;  thus  we  can  expect  a  more  flexible  and  robust 
algorithm.  This  will  be  the  focus  of  the  next  phase  of  this  research. 
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